Lagrangian single particle turbulent statistics through the Hilbert-Huang Transform 
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The Hilbert-Huang transform is applied to analyze single particle Lagrangian velocity data from 
numerical simulations of hydrodynamic turbulence. The velocity trajectory is described in terms 
of a set of intrinsic mode functions, Ci(t), and of their instantaneous frequency, u>i(t). On the 
basis of this decomposition we define the (^-conditioned statistical moments of the Ci modes, named 
g-order Hilbert Spectra (HS). We show that such new quantities have enhanced scaling properties 
as compared to traditional Fourier transform- or correlation-based (Structure Functions) statistical 
indicators, thus providing better insights into the turbulent energy transfer process. We present a 
clear empirical evidence that the energy-like quantity, i.e. the second-order HS, displays a linear 
scaling in time in the inertial range, as expected from dimensional analysis and never observed 
before. We also measure high order moment scaling exponents in a direct way, without resorting the 
Extended Self Similarity (ESS) procedure. This leads to a new estimate of the Lagrangian structure 
functions exponents which are consistent with the multifractal prediction in the Lagrangian frame 
as proposed in [Biferale et al, Phys. Rev. Lett. 93, 064502 (2004)]. 



The statistical description of a tracer trajectory in tur- 
bulent flows still lacks of a sound theoretical and phc- 
nomcnological understanding 0, 0. Presently, no ana- 
lytical results linking the Navier-Stokes equations to the 
statistics of the velocity increments, v(t + r) —v(t), along 
the particle evolution are known. On the ground of di- 
mensional arguments, pure scaling laws are expected for 
time increments larger than the Kolmogorov dissipative 
time, t v , and smaller than the large-scale typical eddy- 
turn over time, 72 . The ratio between the two time scales 
grows with the Reynolds number as Re oc 7l/t j; . Despite 
of the many numerical and experimental attempts [l-Q, 
no clear evidence of scaling properties have been detected 
in the Lagrangian domain even at high Reynolds num- 
bers. Such a fact can be explained cither invoking ultra- 
violet and infrared effects induced by the two cut-offs, t v 
and Tl or by a real pure breaking of scaling invariancc, 
independently of the Reynolds number [H, [j| . Up to now, 
most of the attention has been payed to the so-called 
Lagrangian Structure Functions (LSF), i.e. moments of 
velocity increments: 



the laboratory reference frame. Such relation has been 
well verified in the limit of very small time increments, by 
studying the statistics of flow acceleration [l2[ or by us- 
ing relative scaling properties (l7j . i.e. studying one mo- 
ment versus another one, a procedure known as ESS [l8| . 
On the other hand, no clear evidence of direct scaling 
properties as a function of r has ever been detected (see 
|a |9( for two recent papers discussing this problem) . As 
a result, despite the successful comparisons, using ESS, 
between theoretical predictions for Cl(<?)/Cl(2) and nu- 
merical and experimental Lagrangian measurements (see 
17j). the absence of a clear scaling-range in the time do- 



S g (r) = (\v j (t + T)-v j (t)\ q ), 



(1) 



where for simplicity we have assumed isotropy and 
dropped any possible dependency of the l.h.s on the com- 
ponent of the velocity field. Phenomenological arguments 
based on a 'bridge' relation between Eulerian and La- 
grangian statistics ltj- 16 1 predicts the existence of scal- 
ing properties also in the Lagrangian domain: S q (r) ~ 
r Ctw) for r,, <C r <C 71, with Cl(<?) being related to the 
Eulerian scaling exponents, (e{o.), defining the scaling 
properties of velocity increments between two points in 



main has cast doubts on the one side on the correctness 
and accuracy of the present phenomenological models, 
and on the other side on the fact that SF may not be 
the suitable statistical indicator to study turbulent flows 
in the Lagrangian domain Q . One of the main concern 
regards possible non-local effects induced by either large- 
scales and low-frequencies modes or by small-scales and 
high-frequencies events that may result in sub-leading 
spurious contributions. It is well known for example that 
the temporal evolution of the velocity field along a La- 
grangian trajectory in a turbulent flows is strongly in- 
fluenced by the presence of small-scales vortex filaments 
inducing visible high-frequency oscillations even on the 
single particle velocity signal (sec Fig.Q] and Ref. flij). 
In this paper we want to apply for the first time a rel- 
atively novel technique, called Hilbert-Huang Transform 
(HHT), to analyze multi-scale and multi- frequency sig- 
nals which has revealed to be particularly useful in the 
data analysis of many complex systems HHH. HHT 



has been recently applied to analyze Eulerian turbulent 
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FIG. 1. (a) An example of Lagrangian velocity v(t) with vor- 
tex trapping event from the DNS simulation. The data shows 
the multiscale nature of Lagrangian turbulence with differ- 
ent time scales (structures) superimposed to each other, (b) 
Example of the decomposition of the above trajectory in in- 
trinsic mode functions from empirical mode decomposition. 
Note that the Lagrangian velocity is separated into different 
functions with different time scales. The empirical mode de- 
composition approach reveals the multiscale property of the 
Lagrangian velocity at a local level. 



data [261-1281] , showing an unexpected ability to disentan- 
gle multiscale contributions. The main novelty of HHT 
relies on its frequency-amplitude adaptive nature, being 
based on the decomposition of the original signal on a set 
of quasi- eigenmodes that are not defined a priori [29l. l30j] . 
The idea is to not introduce in the analysis any system- 
atic pre-defined structures as it always happens using 
Fourier-based methodologies (e.g. Fourier decomposition 
or wavelet transforms). 

In this paper, we apply and generalize the HHT 
methodology to extract the hierarchy of Lagrangian scal- 
ing exponent Cl(<?)- The method is applied to the fluid 
trajectory data obtained from Direct Numerical Simu- 
lations (DNS) at Re x = 400 (see Fig.QJ. We present a 
clear empirical evidence of scaling properties in the usual 
sense, as a power of the analyzed frequency, also in the 
Lagrangian domain. We show that the measured Hilbcrt- 
based moments, £ q (u)), display a clear power law on the 
range 0.01 < wt, < 0.1 at least up to the maximum or- 
der allowed to be measured by our statistics, < q < 4. 
The exponents are in good quantitative agreement with 
the one predicted by using the 'bridge relation' based on 
multifractal phenomenology 12(, supporting even more 



the close relationship between Eulerian and Lagrangian 
fluctuations at least for what concerns velocity incre- 
ments in 3D isotropic and homogeneous turbulence. The 
dataset considered here is composed by Lagrangian ve- 
locity trajectories in a homogeneous and isotropic turbu- 
lent flow obtained from a 2048 3 (Re x = 400) DNS sim- 
ulation (more details in [3lj]). We analyze all the avail- 




FIG. 2. Log-log plot of the second order Hilbert Spectrum, 
£2(^0 = E7 = i(|Ci| 2 |cj)t, superposed with the different contri- 
butions from each IMF, {(Ci)' 2 \uj) with i = 1, 6. 

able ~ 2 • 10 5 fluid tracer trajectories, each composed by 
N = 4720 time sampling of Vj{t) (where j = 1,2,3 de- 
notes the three velocity components) saved every 0.1t, ; 
time units. Therefore, we can access time scale from 
0.1 < t/t„ < 236, corresponding to the frequency range 
0.004 < ljt v < 10. The HHT is a procedure composed 
by two steps. The first step is the decomposition of the 
signal into its Intrinsic Mode Functions (IMF) followed 
by the Hilbert transform on such modes. In the first step, 
through a procedure called Empirical Mode Decomposi- 
tion (EMD), we decompose each velocity time series into 
the sum: 



(2) 



where Ci(t) arc the IMFs and r n (t) is a small resid- 
ual, an almost constant function characterized by hav- 
ing at most one extreme along the whole trajectory 
(which will therefore be neglected in the following analy- 
sis) [29|, [33] ■ In eq. ^ n may depend on the trajectory, 
with a maximum value which is linked to its length as 
n ma x = \og 2 (N) ~ 12. Given the actual length of our 
trajectories, with n ~ 6 — 7 we are typically able to re- 
construct the full behaviors (see Fig. [I]). 

To be an IMF, each Ci(t) must satisfy the following 
two conditions: (1) the difference between the number 
of local extrema and the number of zero-crossings must 
be zero or one; (2) the running mean value of the enve- 
lope defined by the local maxima and the envelope de- 
fined by the local minima is zero. Indeed, the IMF is an 
approximation of the so-called mono-component signal, 
whi ch p ossesses a well defined instantaneous frequency 
jiil I32I ]. The physical meaning of such decomposition 
is clear: we want to decompose the original trajectory 
into quasi- eigenmodes with locally homogeneous oscillat- 
ing properties [2^, [33|. In the second step, one performs 
a Hilbert transform for each one of the IMFs, 



(3) 
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where P stands for the Cauchy principal value. This 
allows to retrieve the instantaneous frequency associated 
to each C; via 



Id 

ujAt) — arctan 

n ' 2tt dt 



Cj(t) 

a(t) 



(4) 



|29|. Therefore, we construct the pair of functions 
[Ci(t) , Ui(t)] for all IMF modes, and this concludes the 
standard HHT procedure. Let us stress again the fully 
adaptive nature of the HHT, the IMFs arc not de- 
fined a priori, and they accommodate the oscillatory de- 
gree of the analyzed signal without postulating system- 
atic "structures" [29|, [3(J. The most important conse- 
quence is that the HHT is typically free of sub-harmonics 
2i mill. Here, in order to investigate the amplitude 
of turbulent velocity fluctuations versus their character- 
istic frequency, we define the w-dependent g-order sta- 
tistical moment, C q (u>), by computing the moments of 
each IMF conditioned on those instant of time where the 
corresponding instantaneous frequency has a given value, 
u>i(t) = uj 



£>)=E? =1 (|aHu,) t , 



(5) 



where q > is a real number, and with (. . .) t we denote 
time- and ensemble-averaging over different trajectory re- 
alizations. We dub it Hilbert spectrum (HS) of order q. 
Let us notice that each HS can be seen as a superposition 
of spectra obtained from different IMFs. 

From a dimensional point of view the simplest link 
between the instantaneous frequency uj and the coherence 
time of an eddy r, is the reciprocal relation uj ~ t _1 . 
Therefore, we postulate for the general HS of order q a 
scaling relation of the form 



C q {uj) 



,-a(9) 



(6) 



here, Cl(q) must be compared with the scaling expo- 
nent provided by the LSF (28|. We have validated the 
above scaling relation by using both fractional Brow- 
nian motion with various Hurst number < H < 1 
for mono-fractal processes and a lognormal signal with 
an intermittent parameter [i = 0.15 as an example of 
a multifractal process. For all cases, the scaling expo- 
nents provided by the HHT agree with the ones derived 
by standard SF method and with the theoretical ones 
[281 ] . To begin with, we focus on the case q = 2, that, 
as mentioned, is related to the amplitude of energy fluc- 
tuations as a function of its coherence time or charac- 
teristic frequency. In Fig. [2] we show the second order 
HS, vs uj in log- log, superposed with the contri- 

butions from each different IMF order. As one can see, 
only the whole reconstructed HS shows a good scaling 
behaviour. In order to better compare HS to LSF curves 
we plot them in Fig. [3] in compensated form in such a way 
that the expected behavior in the inertial range would be 




t/t„ l/(/r„), 1/{ut v ) 

FIG. 3. Comparison between the second-order compen- 
sated Lagrangian Structure Function S l 2(r)/(er) vs. t/t v 
(solid line), the compensated Fourier spectrum 
E(f)/e~ 1 f 2 vs. l/{fr v ) (dashed line) and the cor- 
responding Hilbert spectrum C2{uj)€~ 1 uj vs. l/(wr^) 
(•), where t v represents the dissipative time scale of 
the turbulent flow and e the mean energy dissipation 
rate. In the below inset, the logarithmic local slopes for 
dlog&M/dlogr vs. r/rr,, dlog£(/)/dlog/ vs. l/{fr v ) 
and dlog£2(kj)/dloga; vs. 1/(uity;). Note that the expected 
inertial scaling exponents are respectively: 52 (t) ~ r^ 2 ', 
E(f) ~ f-ttiW+V and C 2 (uj) ~ uj- (l(2) , with < L (2) = 1. 




given by a constant, respectively S , 2(r)(er)" 



vs. T 



and 



FIG. 4. Relative contribution of Fourier frequencies in the 
range [f m ,fM] to the ^(t) LSF, as from eq.©. Low 
(IR) frequencies [0, lO -2 ]^ 1 and high (UV) frequencies 
[KF 1 , +cxd]t^ 1 . Vertical lines denote the empirically defined 
inertial range. 

L2(uj)e~ l uj vs. 1/oj. For completeness in the same figure 
also the compensated behavior of the Fourier spectrum, 
E(f)e~ 1 f 2 vs. 1/f, is provided. The first striking dif- 
ference between HS and LSF or Fourier is the enhanced 
scaling property of the new quantity. We also note that 
the shape of LSF curve is consistent with the one in , 
where no plateau was observed in the inertial range. On 
the compensated scale the Fourier spectrum behaves bet- 
ter than the LSF, but the range of scaling is about half 
of that of the Hilbert Spectrum. Such a difference is even 
more evident when the logarithmic local slopes are com- 
pared (see inset of Fig. [3]). A clear inertial scaling range, 
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0.01 < ujT ri < 0.2, corresponding to an interval of time 
scales 5 < t/t v < 100, is observed for the compensated 
£2- The reason why LSF fails in displaying scaling, is 
that it mixes low (infrared, IR)/high (ultraviolet, UV) 
frequency fluctuations to the ones in the inertial-range 
~ [10 _2 ,10 JtjT 1 . This becomes explicit when consid- 
ering the relation, 
and defining 

K f £(T)=S 2 (T)- 



S 2 (r) cc J +o ° E(f)(l - cos27r.fr) df, 



£(/')(!- cos (27r/V))d/', (7) 



which measures the relative contributions to ^(t) from 
the frequency range [f m ,fM]- When such an interval is 
set to [0, lO -2 ]^ 1 we get the low frequency contribu- 
tions, and with [10 _1 , ±oo]t~ 1 the high ones. In figure 
Fig.|4l we show that such spurious non-local contribu- 
tions can be as high as 80%. 




FIG. 5. The Hilbert spectra £ q {u}T v ) for 9 = 1,2,3,4. For 
display clarity, the curves have been vertical shifted by factors 
10" 1 , 10 -2 and 1CF 3 for 9 = 2, 3 and 4. Solid lines comes 
from least square fit in the range is 0.01 < ujt^ < 0.1. The 
inset shows the comparison of the measured local scaling ex- 
ponent (q, oJTrf) — d log C q (uj)/ d log lo with the multifractal 
prediction C,i? F {l)- 

The HS functions C q {uj) have good scaling properties 
also for other q orders. We calculated C q {uj) for the or- 
ders q = 1,2,3,4, and empirically found a good power 
law behavior on the range 0.01 < wt, < 0.1 (resp. 
10 < t/t v < 100), as shown in Fig. [5] This allows to 
extract the scaling exponents directly in the instanta- 
neous frequency space, without resorting to the above 
mentioned ESS procedure. The numerical values for the 
Cl(c?) extracted from the fit in the range 0.01 < ujt v < 0.1 
are reported in the table |U The values of the scaling ex- 
ponents are estimated as the average of the logarithmic 
local slope (L{q 7 u) = dlog£2(w)/dlogw, on the above 
interval and the error bars as the difference between the 
averages taken on only the first or the second half (in 
log scale) of the fitted frequency range. Note that the 



indicated errors are larger than the estimated statisti- 
cal errors. Statistical convergence was here checked by 
performing the same analysis on random subsets with 
1/64 of the total data. First, let us notice the evident 
departure from the dimensional estimate (named K41 
|34|). CF 41 (?) = Second, the measured values are 

in good agreement with the prediction given by the Mul- 
tifractal model, (ff F [Hj]. In order to better appreciate 
the quality of our scaling, we show in the the inset of 
Fig. [5] the logarithmic local slope empirically measured 
with the HHT, £^(9,0;), compensated with the predicted 
value from the multifractal phenomenology, such that a 
plateau around the value 1 is the indication of the exis- 
tence of an intermittent multifractal power law behavior. 





g=l 


g = 2 


9 = 3 


9 = 4 


Cf 41 (9) 


0.5 


1.0 


1.5 


2.0 




0.55 


1 


1.38 


1.71 




0.59 ± 0.06 


1.03 ± 0.03 


1.39 ±0.10 


1.70 ±0.14 



TABLE I. Lagrangian scaling exponents £l (9) for orders 9 = 
1,4 as estimated from dimensional analysis 9/2 (K41), from 
the Multifractal model (MF) [l2T ]. and as obtained here from 
Hilbert Spectra (HS). 



In summary, we have presented a new Hilbert-Huang 
Transform based methodology to capture the intermit- 
tent nature of the turbulent Lagrangian velocity fluctu- 
ations. Our test bench has been a numerical database 
of homogeneous isotropic turbulence at Re\ = 400. The 
first remarkable result is that for the second-order statis- 
tical moment /^(w), an energy-like quantity, we observe 
a clear inertial range versus time defined as r = co^ 1 for 
at least one decade, in the range 0.01 < llit v < 0.2. Such 
clean scaling has never been detected before using more 
standard methods. Second, we extracted the hierarchy of 
scaling exponent Cl(q) for the first time without apply- 
ing ESS. Our measurements provide a solid confirmation 
to the predictions of the multifractal model. The Hilbert 
method we propose in this paper is general and can be 
applied to other systems with multiscale dynamics, e.g., 
Rayleigh-Benard convection (3o| , two dimensional turbu- 
lence (3l,|37|. 
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